function [C,curvelen,xx_int,xx_ext] = get_geometry_moon(ntot,xxc,angle)
% moon
% ntot: number of nodes
% xxc: center of moon
% angle: angle of rotation in radians
rotMat=[cos(angle) -sin(angle);sin(angle) cos(angle)];
tt       = linspace(0,2*pi*(1 - 1/ntot),ntot);
C        = zeros(6,ntot );
C(1,:)   =  cos(tt).*(0.5+cos(tt));
C(2,:)   =  -sin(tt).*(0.5+2*cos(tt));
C(3,:)   =  -cos(tt).*(0.5+2*cos(tt))+2*sin(tt).*sin(tt);
C(4,:)   =   sin(tt);
C(5,:)   =   cos(tt);
C(6,:)   =  -sin(tt);

C([1 4],:)=rotMat*C([1 4],:);
C([2 5],:)=rotMat*C([2 5],:);
C([3 6],:)=rotMat*C([3 6],:);

%approximate length of curve
curvelen=0;
for ii=1:ntot-1
	curvelen=curvelen+sqrt((C(1,ii)-C(1,ii+1))^2+(C(4,ii)-C(4,ii+1))^2);
end

rmin = sqrt(min(C(1,:).^2 + C(4,:).^2));
rmax = sqrt(max(C(1,:).^2 + C(4,:).^2));

%%% Construct interior points.
ttint  = 2*pi*sort(rand(1,3));
xx_int = 0.5*rmin*[cos(ttint); sin(ttint)];

%%% Construct exterior points.
ttext  = 2*pi*sort(rand(1,5));
xx_ext = 1.5*rmax*[cos(ttext); sin(ttext)];

%%% Shift to center at xxc.
C([1,4],:) = C([1,4],:) + xxc*ones(1,size(C,2));
xx_int     = xx_int     + xxc*ones(1,size(xx_int,2));
xx_ext     = xx_ext     + xxc*ones(1,size(xx_ext,2));

return